function[CF1,CF4,CF5,CF6,CF7,CF10] = counterfactual_func(Results,Results_high,CFBaseline,CFBaseline_high,PAR,outputfolder,docounterfactualsims)

PAR.sharenobidders = Results.sharenobidders;

if docounterfactualsims == 1
    FeesBase      = Results;
    FeesBase_high = Results_high;

    % ===================================================== %
    % CF1 : keep cL, cP, cR fixed. Grid of (cB,cS) changes
    % Grid of c
    CF1    = CFBaseline;
    GRIDcB = [-0.35,-0.3,-0.2,-0.1,0,0.1];
    GRIDcS = [-0.1,0,0.1,0.2,0.3,0.35];
    [B,S]  = ndgrid(GRIDcB,GRIDcS);
    Bvec   = B(B>-Inf);
    Svec   = S(S>-Inf);
    x      = length(Bvec);
    shut_bidderentry = 0;
    shut_sellerentry = 0;   

    % Evaluate welfare all fee combinations
    PAR.highvalue = 0;
    for i = 1:x
        Fees = FeesBase;
        Fees.cB = Bvec(i);
        Fees.cS = Svec(i);
        CF1(i)  = Simulate_CF(CFBaseline,PAR,Fees,shut_bidderentry,shut_sellerentry);
    end

    % Do the same for high-value auctions
    CF1_high    = CFBaseline_high;
    PAR.highvalue = 1;
    for i = 1:x
        Fees    = FeesBase_high;
        Fees.cB = Bvec(i);
        Fees.cS = Svec(i);
        CF1_high(i)  = Simulate_CF(CFBaseline_high,PAR,Fees,shut_bidderentry,shut_sellerentry);
    end                       
    
    % ===================================================== %
    % CF4 : specific changes in commissions (to decompose)
    % -- double commission index (cS to 0.204)
    shut_bidderentry = 0;
    shut_sellerentry = 0;
    DoubleCI    = FeesBase;
    DoubleCI.cS = FeesBase.cS*2;
    DoubleCI_high    = FeesBase_high;
    DoubleCI_high.cS = DoubleCI_high.cS*2;
   
    PAR.highvalue = 0;
    CF4  = Simulate_CF(CFBaseline,PAR,DoubleCI,shut_bidderentry,shut_sellerentry);
       
    % ===================================================== %
    % CF5 : Same specific changes in commissions as CF4
    % BUT SHUTTING DOWN SELLER ENTRY
    shut_sellerentry = 1;
    shut_bidderentry = 0;
    PAR.highvalue = 0;
    CF5  = Simulate_CF(CFBaseline,PAR,DoubleCI,shut_bidderentry,shut_sellerentry);
    
    % ===================================================== %
    % CF6 : Same specific changes in commissions as CF4
    % BUT SHUTTING DOWN BIDDER ENTRY
    shut_sellerentry = 0;
    shut_bidderentry = 1;
    PAR.highvalue = 0;
    CF6  = Simulate_CF(CFBaseline,PAR,DoubleCI,shut_bidderentry,shut_sellerentry);
  
    % ===================================================== %
    % CF7 : Same specific changes in commissions as CF4
    % BUT SHUTTING DOWN BIDDER AND SELLER ENTRY
    shut_sellerentry = 1;
    shut_bidderentry = 1;
    PAR.highvalue = 0;
    CF7  = Simulate_CF(CFBaseline,PAR,DoubleCI,shut_bidderentry,shut_sellerentry);   
          
    % ===================================================== %
    % CF9 : Network effects (main sample, positive reserve auctions only)
    % Baseline
    thetas = [Results.mus,Results.ks,Results.shapes];
    thetab = [Results.mub,Results.kb,Results.shapeb];
    PAR.sharenobidders = Results.sharenobidders;
    CDFs  = @(x) ggd_cdf(x,thetas);
    u0vec = unique(sort(arrayfun(@(i) fzero(@(x) CDFs(x)-PAR.v0probs(i),thetas(1)),1:length(PAR.v0probs))));  % simulate seller tastes

    % Baseline values 
    [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,rvec,pipmat,pip0vec,pricemat,highvalmat,reservemat,saleprobmat,price0vec,highval0vec,volumemat,volume0vec] ...
            = piBS_MC(0,0.102,thetab,thetas,PAR);
    T                  = CFBaseline.T*(1-CFBaseline.p_r0);
    M                  = CFBaseline.T*(1-CFBaseline.p_r0)*CFBaseline.lambdastarposr;
    changegrid         = (-100:10:100);
    
    % Value of additional sellers / bidders (network effects) 
    T_grid             = changegrid+CFBaseline.T*(1-CFBaseline.p_r0);
    M_grid             = changegrid+CFBaseline.T*(1-CFBaseline.p_r0)*CFBaseline.lambdastarposr;
    v0bar_grid         = arrayfun(@(i) fsolve(@(x) (T/CFBaseline.CDFu0star)*CDFs(x) - T_grid(i), CFBaseline.u0star),1:length(T_grid));
    Pib_nsgrid         = arrayfun(@(i) Pibidder(Results.eBo_posr,CFBaseline.screeningvalue,v0bar_grid(i),M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    Pib_nsgrid_nov0bar = arrayfun(@(i) Pibidder(Results.eBo_posr,CFBaseline.screeningvalue,CFBaseline.u0star,M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    Pis_nbgrid_marg    = arrayfun(@(i) Piseller(Results.screeningvalue,CFBaseline.CDFu0star,CFBaseline.CDFu0star,Results.eSo+Results.cL,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
    Pis_nbgrid_med     = arrayfun(@(i) Piseller(Results.screeningvalue,CFBaseline.CDFu0star,0.5,Results.eSo+Results.cL,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
   
    % Related change in lambda*, v0* (e.g. out of equilibrium, the other side doesn't adjust yet?)
    lambdastarposr_nsgrid         = arrayfun(@(i) lambdastarfunc(Results.screeningvalue,v0bar_grid(i),Results.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    M_nsgrid                      = T_grid.*lambdastarposr_nsgrid;
    lambdastarposr_nsgrid_nov0bar = arrayfun(@(i) lambdastarfunc(Results.screeningvalue,v0bar_grid(11),Results.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    u0star_nbgrid                 = arrayfun(@(i)  v0starfunc(Results.cL+Results.eSo,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results.screeningvalue,PAR,M_grid(i)/T),1:length(M_grid));   
    CDFu0star_nbgrid              = CDFs(u0star_nbgrid);
    T_nbgrid                      = CDFu0star_nbgrid.*CFBaseline.NS;
    
    % Equilibrium demand elasticity: change lambdastar and vostar when either increasing listing fee or bidder entry fee
    eSgrid = (-2:0.25:2)+(CFBaseline.eSo+CFBaseline.cL);
    eBgrid = (-2:0.25:2)+(CFBaseline.eBo_posr);
    u0star_eSgrid_opt = arrayfun(@(i)  v0starfunc(eSgrid(i),Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results.screeningvalue,PAR),1:length(eSgrid));   
    u0star_eBgrid_opt = arrayfun(@(i)  v0starfunc(Results.cL+Results.eSo,eBgrid(i),Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results.screeningvalue,PAR),1:length(eBgrid));   
    lambdastar_eSgrid_opt = arrayfun(@(i) lambdastarfunc(Results.screeningvalue,u0star_eSgrid_opt(i),Results.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(eSgrid));
    lambdastar_eBgrid_opt = arrayfun(@(i) lambdastarfunc(Results.screeningvalue,u0star_eBgrid_opt(i),eBgrid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(eBgrid));
    
    % ===================================================== %
    % CF9 for high-value too
    % Baseline
    thetas = [Results_high.mus,Results_high.ks,Results_high.shapes];
    thetab = [Results_high.mub,Results_high.kb,Results_high.shapeb];
    PAR.sharenobidders = Results_high.sharenobidders;
    if PAR.generalized==0
        ds    = makedist(PAR.dist,thetas(1),thetas(2));
        CDFs  = @(x) cdf(ds,x);
        v0vec = unique(sort(icdf(PAR.dist,PAR.v0probs,thetas(1),thetas(2))))';         % simulate seller tastes
    else %Generalized gaussian
        CDFs  = @(x) ggd_cdf(x,thetas);
        v0vec = unique(sort(arrayfun(@(i) fzero(@(x) CDFs(x)-PAR.v0probs(i),thetas(1)),1:length(PAR.v0probs))));  % simulate seller tastes
    end
    % Baseline values 
    [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,rvec,pipmat,pip0vec,pricemat,highvalmat,reservemat,saleprobmat,price0vec,highval0vec,volumemat,volume0vec] ...
            = piBS_MC(0,0.09,thetab,thetas,PAR);
    T                  = CFBaseline_high.T*(1-CFBaseline_high.p_r0);
    M                  = CFBaseline_high.T*(1-CFBaseline_high.p_r0)*CFBaseline_high.lambdastarposr;
    changegrid         = (-100:10:100);
    %changegrid         = [-50;-10;0;10;50];
    
    % Value of additional sellers / bidders (network effects) 
    T_grid             = changegrid+CFBaseline_high.T*(1-CFBaseline_high.p_r0);
    M_grid             = changegrid+CFBaseline_high.T*(1-CFBaseline_high.p_r0)*CFBaseline_high.lambdastarposr;
    v0bar_grid         = arrayfun(@(i) fsolve(@(x) (T/CFBaseline_high.CDFu0star)*CDFs(x) - T_grid(i), CFBaseline_high.u0star),1:length(T_grid));
    hPib_nsgrid         = arrayfun(@(i) Pibidder(Results_high.eBo_posr,CFBaseline_high.screeningvalue,v0bar_grid(i),M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    hPib_nsgrid_nov0bar = arrayfun(@(i) Pibidder(Results_high.eBo_posr,CFBaseline_high.screeningvalue,CFBaseline_high.u0star,M/T_grid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    hPis_nbgrid_marg    = arrayfun(@(i) Piseller(Results_high.screeningvalue,CFBaseline_high.CDFu0star,CFBaseline_high.CDFu0star,Results_high.eSo+Results_high.cL,Results_high.eBo_posr,Results_high.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
    hPis_nbgrid_med     = arrayfun(@(i) Piseller(Results_high.screeningvalue,CFBaseline_high.CDFu0star,0.5,Results_high.eSo+Results_high.cL,Results_high.eBo_posr,Results_high.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,M_grid(i)/T),1:length(M_grid));
   
    % Related change in lambda*, v0* (e.g. out of equilibrium, the other side doesn't adjust yet?)
    hlambdastarposr_nsgrid         = arrayfun(@(i) lambdastarfunc(Results_high.screeningvalue,v0bar_grid(i),Results_high.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    hM_nsgrid                      = T_grid.*hlambdastarposr_nsgrid;
    hlambdastarposr_nsgrid_nov0bar = arrayfun(@(i) lambdastarfunc(Results_high.screeningvalue,v0bar_grid(11),Results_high.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(T_grid));
    u0star_nbgrid                 = arrayfun(@(i)  v0starfunc(Results_high.cL+Results_high.eSo,Results_high.eBo_posr,Results_high.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results_high.screeningvalue,PAR,M_grid(i)/T),1:length(M_grid));   
    CDFu0star_nbgrid              = CDFs(u0star_nbgrid);
    hT_nbgrid                      = CDFu0star_nbgrid.*CFBaseline_high.NS;
    
    % Equilibrium demand elasticity: change lambdastar and vostar when either increasing listing fee or bidder entry fee
    eSgrid = (-2:0.25:2)+(CFBaseline_high.eSo+CFBaseline_high.cL);
    eBgrid = (-2:0.25:2)+(CFBaseline_high.eBo_posr);
    u0star_eSgrid_opt = arrayfun(@(i)  v0starfunc(eSgrid(i),Results_high.eBo_posr,Results_high.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results_high.screeningvalue,PAR),1:length(eSgrid));   
    u0star_eBgrid_opt = arrayfun(@(i)  v0starfunc(Results_high.cL+Results_high.eSo,eBgrid(i),Results_high.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,Results_high.screeningvalue,PAR),1:length(eBgrid));   
    hlambdastar_eSgrid_opt = arrayfun(@(i) lambdastarfunc(Results_high.screeningvalue,u0star_eSgrid_opt(i),Results_high.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(eSgrid));
    hlambdastar_eBgrid_opt = arrayfun(@(i) lambdastarfunc(Results_high.screeningvalue,u0star_eBgrid_opt(i),eBgrid(i),pibmat,pib0vec,pdfmat,v0vec,0,PAR), 1:length(eBgrid));

    
    % ===================================================== %
   
% CF10 Effect of increasing fees across distribution of V0 --> lemon story
    % Grid of v0|v0 r>0 (main sample)
    quantilegrid = (0.1:0.1:0.9);
    PAR.highvalue=0;
    
    % Baseline values
    thetas = [Results.mus,Results.ks,Results.shapes];
    thetab = [Results.mub,Results.kb,Results.shapeb];
    PAR.sharenobidders = Results.sharenobidders;
    CDFs  = @(x) ggd_cdf(x,thetas);
    [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec] = piBS_MC(0,0.102,thetab,thetas,PAR);
  
    % Baseline surplus 
    Pis_v0grid                  = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CFBaseline.CDFu0star,quantilegrid(i),Results.eSo+Results.cL,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,CFBaseline.lambdastarposr),1:length(quantilegrid));

    % Increase cL (+1)
    Pis_v0grid_plus1cL_noentry  = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CFBaseline.CDFu0star,quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,CFBaseline.lambdastarposr),1:length(quantilegrid));
    v0star_plus1cL              = v0starfunc(Results.cL+Results.eSo+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,CFBaseline.screeningvalue,PAR);
    Pis_v0grid_plus1cL          = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CDFs(v0star_plus1cL),quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));

    % how many are pushed out?
    (CFBaseline.NS*CDFs(v0star_plus1cL)-CFBaseline.NS*CFBaseline.CDFu0star)/ (CFBaseline.NS*CFBaseline.CDFu0star)

    % how does this affect number of bidders / listing?
    m1lambdastar                  = lambdastarfunc(CFBaseline.screeningvalue,v0star_plus1cL,Results.eBo_posr,pibmat,pib0vec,pdfmat,v0vec,0,PAR);
    (CFBaseline.lambdastarposr-m1lambdastar)/CFBaseline.lambdastarposr
    



% --- m2: More seller heterogeneity
    % Variance distribution: 
    varGG = @(thetas) ((thetas(2)^2)/(thetas(3)^2))* exp(thetas(3)^2)*( exp(thetas(3)^2)-1);
    oldVariance=varGG(thetas);
    
    thetas = [Results.mus*1.04,Results.ks*1.4,Results.shapes];%to increase sigma by 40%
    newVariance=varGG(thetas); % New variance
    newVariance/oldVariance; %Increases variance by this much
    thetab = [Results.mub,Results.kb,Results.shapeb];
    PAR.sharenobidders = Results.sharenobidders;
    if PAR.generalized==0
        ds    = makedist(PAR.dist,thetas(1),thetas(2));
        CDFs  = @(x) cdf(ds,x);
    else %Generalized gaussian
        CDFs  = @(x) ggd_cdf(x,thetas);
    end
    
    % Baseline values
    [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec] = piBS_MC(0,0.102,thetab,thetas,PAR);   

    % Baseline surplus 
    m2v0star                      = v0starfunc(Results.cL+Results.eSo,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,CFBaseline.screeningvalue,PAR);
    m2Pis_v0grid                  = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CDFs(m2v0star),quantilegrid(i),Results.eSo+Results.cL,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));
    
    % Increase cL (+1)
    m2v0star_plus1cL              = v0starfunc(Results.cL+Results.eSo+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,CFBaseline.screeningvalue,PAR);
    m2Pis_v0grid_plus1cL          = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CDFs(m2v0star_plus1cL),quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));
      
    % how many are pushed out with added heterogeneity? Calibrated to be
    % the same share as line 189 for comparable statements (~5.5 percent)
    (CFBaseline.NS*CDFs(m2v0star_plus1cL)-CFBaseline.NS*CDFs(m2v0star))/ (CFBaseline.NS*CDFs(m2v0star));


% --- m3: Use the +1 pound to subsidize bidder entry
    thetas = [Results.mus,Results.ks,Results.shapes];
    thetab = [Results.mub,Results.kb,Results.shapeb];
    PAR.sharenobidders = Results.sharenobidders;
    if PAR.generalized==0
        ds    = makedist(PAR.dist,thetas(1),thetas(2));
        CDFs  = @(x) cdf(ds,x);
    else %Generalized gaussian
        CDFs  = @(x) ggd_cdf(x,thetas);
    end
    
    % Baseline values
    [pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec]  = piBS_MC(0,0.102,thetab,thetas,PAR);
          
    % Increase cL (+1), subsidy giving extra revenue back to
    % entrants -- value chosen to keep the marginal seller indifferent
    subsidy = 0.18;
    m3Pis_v0grid_plus1cL_noentry  = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CFBaseline.CDFu0star,quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,CFBaseline.lambdastarposr),1:length(quantilegrid));
    m3v0star_plus1cL              = v0starfunc(Results.cL+Results.eSo+1,Results.eBo_posr-subsidy,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,CFBaseline.screeningvalue,PAR);
    m3lambdastar                  = lambdastarfunc(CFBaseline.screeningvalue,m3v0star_plus1cL,Results.eBo_posr-subsidy,pibmat,pib0vec,pdfmat,v0vec,0,PAR);
    m3Pis_v0grid_plus1cL          = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CDFs(m3v0star_plus1cL),quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr-subsidy,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));
    
    % Marginal seller indifferent condition: must be ~0
    m3v0star_plus1cL-CFBaseline.u0star

    % Effect on platform revenues and other groups
    Fees = CFBaseline;
    Fees.cB=0;
    Fees.cS=0.102;
    Fees.cL=3.1;
    Fees.cP=-subsidy;
    CF10 = Simulate_CF(CFBaseline,PAR,Fees,0,0);
    (CF10.P-CFBaseline.P)/CFBaseline.P;
    (CF10.Volume-CFBaseline.Volume)/CFBaseline.Volume;
    (CF10.EWB-CFBaseline.EWB)/CFBaseline.EWB;
    (CF10.EWB_posr-CFBaseline.EWB_posr)/CFBaseline.EWB_posr;
    (CF10.EWB_nor-CFBaseline.EWB_nor)/CFBaseline.EWB_nor; % Even sellers r=0 better off
    (CF10.ES_posr-CFBaseline.ES_posr)/CFBaseline.ES_posr;
    (CF10.ES_nor-CFBaseline.ES_nor)/CFBaseline.ES_nor;
    (CF10.T-CFBaseline.T)/CFBaseline.T;
    
    
% --- Budget neutral bidder subsidy
    Fees = CFBaseline;
    Fees.cB=0;
    Fees.cS=0.102;
    Fees.cL=3.1;
    BNsubsidy=0.23;
    Fees.cP=-BNsubsidy;
    CF10b = Simulate_CF(CFBaseline,PAR,Fees,0,0);
    (CF10b.P-CFBaseline.P)/CFBaseline.P %BUDGET NEUTRALITY CONDITION: this should be ~0

    BNm3Pis_v0grid_plus1cL_noentry  = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CFBaseline.CDFu0star,quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR,CFBaseline.lambdastarposr),1:length(quantilegrid));
    BNm3v0star_plus1cL              = v0starfunc(Results.cL+Results.eSo+1,Results.eBo_posr-BNsubsidy,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,CFBaseline.screeningvalue,PAR);
    BNm3lambdastar                  = lambdastarfunc(CFBaseline.screeningvalue,m3v0star_plus1cL,Results.eBo_posr-BNsubsidy,pibmat,pib0vec,pdfmat,v0vec,0,PAR);
    BNm3Pis_v0grid_plus1cL          = arrayfun(@(i) Piseller(CFBaseline.screeningvalue,CDFs(m3v0star_plus1cL),quantilegrid(i),Results.eSo+Results.cL+1,Results.eBo_posr-BNsubsidy,Results.cR,thetas,pibmat,pib0vec,pdfmat,pismat,pis0mat,v0vec,0,PAR),1:length(quantilegrid));
        

   % Save results
    save(strcat(outputfolder,'Counterfactual_results.mat'))
else
    load(strcat(outputfolder,'Counterfactual_results.mat'))
end

    %% Some useful baseline statistics
    X = CFBaseline;
       
    % Baseline main:
    X.P      = CFBaseline.P;
    X.Volume = CFBaseline.Volume;
    X.EWB    = CFBaseline.EWB; 
    X.ES     = CFBaseline.ES; 
    X.TOT    = X.EWB+X.ES;
    X.T                 = CFBaseline.T;
    X.ESperseller       = X.ES/X.T;
    X.EWBperbidder      = X.EWB/X.T;
    X.winbid = CFBaseline.winbid;
    X.saleprob = CFBaseline.saleprob;

    % Baseline high:
    X.P_high      = CFBaseline_high.P;
    X.Volume_high = CFBaseline_high.Volume;
    X.EWB_high    = CFBaseline_high.EWB; 
    X.ES_high     = CFBaseline_high.ES; 
    X.TOT_high    = X.EWB_high+X.ES_high;
    X.T_high            = CFBaseline_high.T;
    X.ESperseller_high  = X.ES_high/X.T_high;
    X.EWBperbidder_high = X.EWB_high/X.T_high;
    X.winbid_high = CFBaseline_high.winbid;
    X.saleprob_high = CFBaseline_high.saleprob;
    X.CDFu0star_high = CFBaseline_high.CDFu0star;
    X.lambdastarposr_high = CFBaseline_high.lambdastarposr;

    % Share high
    X.sharehighvalueT = CFBaseline_high.T/(CFBaseline_high.T+CFBaseline.T);
    X.sharehighvalueM = CFBaseline_high.M/(CFBaseline_high.M+CFBaseline.M);
    X.sharehighvolume = CFBaseline_high.Volume/(CFBaseline_high.Volume + CFBaseline.Volume);
    X.sharehighprofit = CFBaseline_high.P/(CFBaseline_high.P+CFBaseline.P);
 
    
    %% ===================================================== %
    % CF1 : keep cL, cP, cR fixed. Grid of (cB,cS) changes
    GRIDcB = [-0.35,-0.3,-0.2,-0.1,0,0.1];
    GRIDcS = [-0.1,0,0.1,0.2,0.3,0.35];
    [B,S]  = ndgrid(GRIDcB,GRIDcS);
    Bvec   = B(B>-Inf);
    Svec   = S(S>-Inf);
    x      = length(Bvec);
    
    % Create a grid and a finer fee grid
    fineGRIDcS      = linspace(min(GRIDcS),max(GRIDcS),41);
    fineGRIDcB      = linspace(min(GRIDcB),max(GRIDcB),41);
    [Sq,Bq]         = meshgrid(Svec,Bvec);                               % To create the 2D matrix of coarse cS,cB coordinates
    [Sqfine,Bqfine] = meshgrid(fineGRIDcS,fineGRIDcB);                   % Fine mesh grid
    
    % Get Platform Revenue, User Welfare, main / high / combined [assuming
    % share of high-value auctions remains the same in CF]
    CF1_PlatRev  = NaN(x,1);
    CF1_SelSurp  = NaN(x,1);
    CF1_BidSurp  = NaN(x,1);
    CF1_Volume   = NaN(x,1);
    CF1_Users    = NaN(x,1);
    CF1_TotalRev = NaN(x,1);
    CF1_UserSurp = NaN(x,1);
    CF1_lambdastarposrMAIN = NaN(x,1);
    CF1_CDFu0starMAIN = NaN(x,1);
    CF1_VolumeMAIN = NaN(x,1);
    CF1_PlatRevMAIN = NaN(x,1);

    for i=1:x %add CF1_high to all results
        CF1(i).CDFu0star ;
        if CF1(i).CDFu0star < 0.999
            CF1_PlatRev(i)  = CF1(i).P + CF1_high(i).P;
            CF1_SelSurp(i)  = CF1(i).ES + CF1_high(i).ES;
            CF1_BidSurp(i)  = CF1(i).EWB + CF1_high(i).EWB;
            CF1_Volume(i)   = CF1(i).Volume + CF1_high(i).Volume;
            CF1_Users(i)    = CF1(i).T+CF1(i).M + CF1_high(i).T + CF1_high(i).M;
            CF1_TotalRev(i) = CF1(i).Total + CF1_high(i).Total;
            CF1_UserSurp(i) = CF1(i).EWB+ CF1(i).ES + CF1_high(i).EWB + CF1_high(i).ES;
            
            CF1_CDFu0starMAIN(i) = CF1(i).CDFu0star ;
            CF1_VolumeMAIN(i) = CF1(i).Volume;
            CF1_PlatRevMAIN(i) = CF1(i).P';
            CF1_lambdastarposrMAIN(i) = CF1(i).lambdastarposr;
        else
            CF1_PlatRev(i)  = NaN;
            CF1_SelSurp(i)  = NaN;
            CF1_BidSurp(i)  = NaN;
            CF1_Volume(i)   = NaN;
            CF1_Users(i)    = NaN;
            CF1_TotalRev(i) = NaN;
            CF1_UserSurp(i) = NaN;
            
            CF1_CDFu0starMAIN(i) = NaN ;
            CF1_lambdastarposrMAIN(i)= NaN;
            CF1_PlatRevMAIN(i) = NaN;
            CF1_VolumeMAIN(i) = NaN;
        end
    end
    
      
% ----- Contour plots to illustrate Commission Index ---- %
    % Main only: for clean analysis
        
    % Platform Revenue & commission index
    F          = scatteredInterpolant(Svec,Bvec,CF1_PlatRevMAIN);
    MAINPlatRev_cScB     = F(Sq,Bq);                                                 
    MAINPlatRev_cScBfine = griddata(Sq,Bq,MAINPlatRev_cScB,Sqfine,Bqfine,'linear');  
    MAINPlatRev_cScBfineshare = MAINPlatRev_cScBfine / (CFBaseline.P);
    
     % CDFu0star & commission index
    F          = scatteredInterpolant(Svec,Bvec,CF1_VolumeMAIN);
    MAINVolume_cScB     = F(Sq,Bq);                                                 
    MAINVolume_cScBfine = griddata(Sq,Bq,MAINVolume_cScB,Sqfine,Bqfine,'linear');  
    MAINVolume_cScBfineshare = MAINVolume_cScBfine / (CFBaseline.Volume);
    
    % Generate commission index
    ComIdx  = (Bqfine+Sqfine)./(1+Bqfine);
    ComIdxshare  = ComIdx  ./ ((CFBaseline.cB+CFBaseline.cS)./(1+CFBaseline.cB));
    
 
    % THIS ONE Fig 2a
    % Plot for Volume
    figure
    contourf(Sqfine,Bqfine,MAINVolume_cScBfineshare,[0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8],'ShowText','on')
    xlabel('Seller commission')
    ylabel('Buyer commission')
    ax = gca;
    ax.FontSize=14;
    ax.TitleFontSizeMultiplier = 1.2;
    axis square
    colormap(gray(15))
    brighten(0.5)
    hold on
    contour(Sqfine,Bqfine,ComIdxshare,[0.5 0.5],'LineWidth',7,'LineStyle','--','color','black','ShowText','on')
    contour(Sqfine,Bqfine,ComIdxshare,[1 1],'LineWidth',7,'LineStyle', '--','color','black','ShowText','on')
    contour(Sqfine,Bqfine,ComIdxshare,[2 2],'LineWidth',7,'LineStyle', '--', 'color','black','ShowText','on')
    line([-0.1 0.35],[0 0], 'LineWidth',1,'Color','black')
    line([0 0],[-0.35 0.1], 'LineWidth',1,'Color','black')
    plot(0.102,0,'wpentagram','MarkerSize',20,'MarkerFaceColor','white')
    xlim([-0.1 0.3])
    ylim([-0.3 0.1])
    legend('Simulated Volume levels','Commission index levels','Location','SouthWest')
    print([outputfolder,'Figure2a'],'-dpng')

    % Plot for Platform Revenue - Fig 2b
    figure
    contourf(Sqfine,Bqfine,MAINPlatRev_cScBfineshare,[0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8],'ShowText','on')
    xlabel('Seller commission')
    ylabel('Buyer commission')
    ax = gca;
    ax.FontSize=14;
    ax.TitleFontSizeMultiplier = 1.2;
    axis square
    colormap(gray(15))
    brighten(0.5)
    hold on
    contour(Sqfine,Bqfine,ComIdxshare,[1 1],'LineWidth',7,'LineStyle','--','color','black','ShowText','on')
    contour(Sqfine,Bqfine,ComIdxshare,[2 2],'LineWidth',7, 'LineStyle', '--', 'color','black','ShowText','on')
    contour(Sqfine,Bqfine,ComIdxshare,[0.4 0.4],'LineWidth',7,'LineStyle','--','color','black','ShowText','on')
    line([-0.1 0.35],[0 0], 'LineWidth',1,'Color','black')
    line([0 0],[-0.35 0.1], 'LineWidth',1,'Color','black')
    plot(0.102,0,'wpentagram','MarkerSize',20, 'MarkerFaceColor','white')
    xlim([-0.1 0.3])
    ylim([-0.3 0.1])
    legend('Simulated Platform Revenue levels','Commission index levels','Location','SouthWest')
    print([outputfolder,'Figure2b'],'-dpng')

    
    % ---- Other contour plots ---- %   
    % Pre-calculate values for contour plots
    % Revenue
    F          = scatteredInterpolant(Svec,Bvec,CF1_PlatRev);
    P_cScB     = F(Sq,Bq);                                                % To evaluate one-D outcome (PcP0) on 2D grid
    P_cScBfine = griddata(Sq,Bq,P_cScB,Sqfine,Bqfine,'linear');           % Evaluate function (direct revenues conditional) on grid
    P_cScBfineshare = P_cScBfine ./ (CFBaseline.P+CFBaseline_high.P);     % Share of current

    % Volume for contour line at baseline level
    F = scatteredInterpolant(Svec,Bvec,CF1_Volume);
    V_cScB = F(Sq,Bq);
    V_cScBfine = griddata(Sq,Bq,V_cScB,Sqfine,Bqfine,'linear');    
    V_cScBfineshare = V_cScBfine ./ (CFBaseline.Volume+CFBaseline_high.Volume);

    % Fig 4
    % Platform Revenue with volume constraint (including high-value sample)
    figure
    contourf(Sqfine,Bqfine,P_cScBfineshare,[0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8],'ShowText','on')
    xlabel('Seller commission')
    ylabel('Buyer commission')
    ax = gca;
    ax.FontSize=14;
    ax.TitleFontSizeMultiplier = 1.2;
    axis square
    colormap(gray(15))
    brighten(0.5)
    hold on
    contour(Sqfine,Bqfine,V_cScBfineshare,[1 1],'LineWidth',7,'color','black','LineStyle', '--','ShowText','on')
    fill([-0.1,-0.1,0.35,0.35],[-0.003,0.003,0.003,-0.003],'black','LineStyle','none')
    alpha(0.2)
    line([-0.1 0.35],[0 0], 'LineWidth',1,'Color','black')
    line([0 0],[-0.25 0.1], 'LineWidth',1,'Color','black')
    fill([CFBaseline.cS,CFBaseline.cS,CFBaseline_high.cS,CFBaseline_high.cS],[-0.35,0.1,0.1,-0.35],'black','LineStyle','none')
    plot(0.098,0,'wpentagram','MarkerSize',20, 'MarkerFaceColor','white')
    alpha(0.2)
    xlim([-0.1 0.35])
    ylim([-0.22 0.1])
    legend('Simulated Platform Revenue levels','Volume constraint','Location','SouthWest')
    print([outputfolder,'Figure4'],'-dpng')
     
         
   
    %% ===================================================== %
    % CF9 : Simulated Network Effects
    % main sample and high-end samples, positive reserve auctions
 
    % Table 6

    % For table: indirect network effect seller entry on bidders
    p=polyfit(changegrid, Pib_nsgrid,4);
    a=polyval(p,changegrid);
    a=a-a(11); % 100% selection
    b= lambdastarposr_nsgrid-lambdastarposr_nsgrid(11); % Eqm response lambdastar
    c=     M_nsgrid   - M_nsgrid(11); % Eqm response M
    d=(0.9.*polyval(p,changegrid)+0.1.*Pib_nsgrid_nov0bar); % 90% selection
    d=d-d(11);
    e=Pib_nsgrid_nov0bar-Pib_nsgrid_nov0bar(11); % No selection
  
    % For table: indirect network effect bidder entry on sellers
    f=Pis_nbgrid_marg-Pis_nbgrid_marg(11);
    g= u0star_nbgrid-u0star_nbgrid(11); % Eqm v0star
    h=    T_nbgrid   - T_nbgrid(11); % Eqm T
    i=Pis_nbgrid_med-Pis_nbgrid_med(11);
    
    % High end
    % Indirect network effect seller entry on bidders
    p=polyfit(changegrid, hPib_nsgrid,4);
    ha=polyval(p,changegrid);
    ha=ha-ha(11); % 100% selection
    hd=(0.9.*polyval(p,changegrid)+0.1.*hPib_nsgrid_nov0bar); % 90% selection
    hd=hd-hd(11);
    he=hPib_nsgrid_nov0bar-hPib_nsgrid_nov0bar(11); % No selection
  
    % For table: indirect network effect bidder entry on sellers
    hf=hPis_nbgrid_marg-hPis_nbgrid_marg(11);
    hi=hPis_nbgrid_med-hPis_nbgrid_med(11);

    
     % Smaller table 
     %- effect sellers on bidders
    indx = [6,10,12,16];
    sellersonbidders = [a(indx);e(indx)];
     %- effect bidders on sellers
    biddersonsellers = [f(indx);i(indx)];
    % high end
    indx = [6,16];
    hsellersonbidders = [ha(indx);he(indx)];
     %- effect bidders on sellers
    hbiddersonsellers = [hf(indx);hi(indx)];
   

    matrix2latex(round([[sellersonbidders;biddersonsellers],[hsellersonbidders;hbiddersonsellers]],3),...
    [outputfolder,'Table6.tex'],...
    'ColumnLabels',[{'-50'},{'-10'},{'+10'},{'+50'},{'-50h'},{'+50h'}],...
    'RowLabels',[{'Selection'},{'No Selection'},{'Marginal seller'},{'Median seller'}])
 
 
    
    %% Lemon story (plotted for inframarginal sellers)
    % Figure 3a
    f=figure
    xlabel('Decile seller value distribution')
    ylabel('Change seller surplus')
    f.Position = [10 10 550 550]; 
    ax = gca;
    ax.FontSize=14;
    ax.TitleFontSizeMultiplier = 1.2;
    axis square
    hold on
    line(quantilegrid,Pis_v0grid_plus1cL_noentry-Pis_v0grid,'LineWidth',1,'color','black','Marker','o','MarkerSize',10)
    line(quantilegrid,Pis_v0grid_plus1cL-Pis_v0grid,'LineWidth',2,'color','black','Marker','s','MarkerSize',10)
    line(quantilegrid,m2Pis_v0grid_plus1cL-m2Pis_v0grid,'LineWidth',2,'color','black','Marker','x','MarkerSize',10)
    line([0.09 0.91],[0 0],'LineStyle','--','LineWidth',0.5,'color','black')
    xlim([0.09, CFBaseline.CDFu0star])
    ylim([-1.1 -0.4])
    legend('No entry','Two-sided entry','Added seller heterogeneity','Location','NorthEast')
    print([outputfolder,'Figure3a'],'-dpng')

    % Figure 3b
    % What about (budget-neutral) changes in entry fees
    f=figure
    xlabel('Decile seller value distribution')
    ylabel('Change seller surplus')
    ax = gca;
    ax.FontSize=14;
    ax.TitleFontSizeMultiplier = 1.2;
    f.Position = [10 10 550 550]; 
    axis square
    hold on
    line(quantilegrid,Pis_v0grid_plus1cL_noentry-Pis_v0grid,'LineWidth',1,'color','black','Marker','o','MarkerSize',10)
    line(quantilegrid,m3Pis_v0grid_plus1cL-Pis_v0grid,'LineWidth',2,'color','black','Marker','s','MarkerSize',10)%marg seller indifferent
    line(quantilegrid,BNm3Pis_v0grid_plus1cL-Pis_v0grid,'LineWidth',2,'color','black','Marker','o','MarkerSize',10)%Budget neutral
    line([0.09 0.9],[0 0],'LineStyle','--','LineWidth',0.5,'color','black')
    xlim([0.09 CFBaseline.CDFu0star])
    ylim([-1.1 2.5])
    legend('No entry','Two-sided entry + marg seller-indifferent bidder subsidy','Two-sided entry + budget-neutral bidder subsidy','Location','NorthEast')
    print([outputfolder,'Figure3b'],'-dpng')
    
    
    
    %% Damages effects referenced in text
    % Part I: double commission index: cS*2
    % CF4: two-sided entry, CF5: shutting down seller entry, CF6: shutting
    % down bidder entry, CF7: shutting down both
    Total_damage   = -1*[(CF7.EWB_posr+CF7.ES_posr)-(CFBaseline.EWB_posr+CFBaseline.ES_posr); %shutting down entry
                     (CF5.EWB_posr+CF5.ES_posr)-(CFBaseline.EWB_posr+CFBaseline.ES_posr);  %shutting down seller entry
                     (CF4.EWB_posr+CF4.ES_posr)-(CFBaseline.EWB_posr+CFBaseline.ES_posr)];  %two-sided entry
    Hammer_change  = [(CF7.winbid_posr*CF7.saleprob_posr - CFBaseline.winbid_posr*CFBaseline.saleprob_posr)/(CFBaseline.winbid_posr*CFBaseline.saleprob_posr);
                      (CF5.winbid_posr*CF5.saleprob_posr - CFBaseline.winbid_posr*CFBaseline.saleprob_posr)/(CFBaseline.winbid_posr*CFBaseline.saleprob_posr);
                      (CF4.winbid_posr*CF4.saleprob_posr - CFBaseline.winbid_posr*CFBaseline.saleprob_posr)/(CFBaseline.winbid_posr*CFBaseline.saleprob_posr)];
    Tot_post_ham   = [CF7.T*(1-CF7.p_r0)*CF7.winbid_posr*CF7.saleprob_posr;
                                    CF5.T*(1-CF5.p_r0)*CF5.winbid_posr*CF5.saleprob_posr;
                                    CF4.T*(1-CF4.p_r0)*CF4.winbid_posr*CF4.saleprob_posr];
    Seller_damage  = -1*[(CF7.ES_posr-CFBaseline.ES_posr);
                      (CF5.ES_posr-CFBaseline.ES_posr);
                      (CF4.ES_posr-CFBaseline.ES_posr)];
    Incidence_sel  = Seller_damage./Total_damage        ;      
    Damage_per_wb_postham  = [Total_damage-Seller_damage]./Tot_post_ham  ;    
    Damage_per_sel_postham = [Seller_damage]./Tot_post_ham   ;

    matrix2latex([round([Total_damage/100000,Incidence_sel,Hammer_change,Damage_per_wb_postham,Damage_per_sel_postham]*100,1)]',...
    [outputfolder,'Damages_table.tex'],...
    'RowLabels',[{'Total damage'},{'Incidence sellers'},{'Change hammer price'},{'E buyer damage (share hammer)'},{'E seller damage (share hammer)'}],...
    'ColumnLabels',[{'No entry'},{'No seller entry'},{'Full two-sided entry'}])

    % Stat in paper: increase in seller damage: 
    (Seller_damage(3)-Seller_damage(1))/(Seller_damage(1))

    
end

